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Abstract 

Modeling  of  multi-component  gas  transport  (H2,  N2,  H20)  in  a  solid  oxide  fuel  cell  (SOFC)  anode  has  been  carried  out  using  a  recently 
developed  two-dimensional  (2D)  lattice  Boltzmann  method  (LBM)  model.  The  LBM  can  simulate  gas  diffusion  in  complex  porous  media  between 
species  having  different  molecular  weights.  The  porous  anode  structure  through  which  gas  transport  occurs  is  obtained  from  micrograph  images 
of  an  existing  SOFC,  converted  to  digital  form  and  used  as  an  input  for  the  LBM  model.  Effect  of  medium  geometry,  mainly  the  porosity,  and 
dimensionless  flux  on  the  fuel  (H2)  delivery  to  the  active  sites  and  product  (H20)  removal  is  examined.  The  relationship  between  mass  transfer  and 
concentration  polarization  in  the  anode  has  been  clarified.  Predicted  concentration  polarization  plots  using  the  LBM  model  are  validated  against 
prior  studies  and  show  good  agreement.  Because  of  its  ability  to  incorporate  detailed  geometry  information,  the  LBM  model  can  be  used  for  design 
and  optimization  of  SOFC  electrodes  without  empirical  modification  of  diffusion  coefficients  using  medium  porosity  and  tortuosity.  As  more 
advanced  three-dimensional  anode  imaging  techniques  are  developed,  the  LBM  model  can  be  adapted  to  model  current  flow  through  the  anode 
material  in  addition  to  mass  transport  through  the  anode  pores. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

In  a  typical  solid  oxide  fuel  cell  (SOFC)  [Fig.  1(a)],  an  elec¬ 
trochemical  reaction  described  by  Eq.  (1)  occurs  at  the  anode 
triple  phase  boundary  (TPB).  The  anode  TPB  is  the  interface 
between  the  porous  anode  material  (nickel  yttria- stabilized  zir- 
conia,  Ni  +  YSZ),  electrolyte  (YSZ)  and  the  gas  species  (H2  and 
H2O)  present  in  the  pores.  A  hydrogen  (H2)  molecule  diffusing 
through  the  pores  in  the  anode  from  the  fuel  or  gas  side  at  x  =  0  to 
the  TPB  at  x-L  combines  with  a  negatively  charged  oxygen  ion 
(O2-)  traveling  through  the  solid  electrolyte  from  the  cathode 
side.  H2  forms  a  molecule  of  H2O  during  this  oxidation  reaction 
and  releases  a  pair  of  electrons 

H2+02“^  H2O  +  2e“  (1) 

l/202  +  2e“  — >  O2-  (2) 
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The  H2O  molecule  formed  at  the  TPB  then  diffuses  back 
through  the  pores  in  the  anode  from  a  =  L  towards  the  gas  channel 
at  x  =  0.  The  electrons  generated  at  the  anode  TPB  are  con¬ 
ducted  through  the  anode  material  to  the  current  collector,  and 
flow  through  an  external  circuit.  The  electrons  finally  enter  the 
cathode  and  are  conducted  through  the  cathode  material  to  the 
cathode  TPB.  Here,  the  electrons  combine  with  oxygen  that  has 
diffused  from  the  air  through  pores  in  the  cathode  to  create  O2- 
ions  as  indicated  in  Eq.  (2).  These  ions  are  conducted  through 
the  electrolyte  and  the  cycle  is  completed.  Assuming  that  each 
hydrogen  molecule  is  consumed  and  forms  a  water  molecule  at 
the  anode  TPB,  the  mole  flux  of  hydrogen  is  equal  and  opposite 
to  the  mole  flux  of  water  vapor.  This  case  is  often  referred  to 
as  equimolar  counter-diffusion.  Actual  diffusion  in  the  porous 
anode  occurs  in  the  presence  of  a  third  species,  nitrogen  (N2). 
The  N2  takes  no  active  part  in  the  chemical  reaction;  however 
the  presence  of  N2  affects  transport  of  other  species  within  the 
fuel  cell. 

The  goal  of  this  work  is  to  model  the  diffusive  flow  of  three 
species  (H2,  H2O  and  N2)  through  the  porous  anode  structure 
using  a  two-dimensional  (2D)  representation,  so  that  transport 
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Nomenclature 

Ct  total  molar  concentration  (mol  m-3) 

Du  binary  diffusivity  between  H2  and  H2O  (m2  s-1) 

£>23  binary  diffusivity  between  H2O  and  N2  (m2  s-1) 

£>3i  binary  diffusivity  between  N2  and  H2  (m2  s-1) 
e  base  velocity  in  the  lattice  Boltzmann  model 

E  actual  voltage  developed  by  the  fuel  cell  (V) 

Eo  maximum  voltage  (V) 

/  particle  velocity  distribution  function 

F  Faradays  constant  (96,485  A  s  mol- 1 ) 

H  height  of  the  model  anode  (m) 

i  current  density  (A  m-2) 

t  non-dimensional  current  density 

J  mole  flux  (mol  m-2  s- 1 ) 

J*  dimensionless  flux  =  ££/Ct£>23 

L  total  thickness  of  the  anode  (m) 

N  number  of  processors  used  in  a  parallel  computa¬ 
tion 

P  partial  pressure 

R  universal  gas  constant  (8.3145  J  mol-1  K-1) 

t  time  (s) 

v  distance  in  the  anode,  measured  from  the  gas  side 

(m) 

x  position  of  a  lattice  node 

v*  distance  normalized  by  the  total  length 

X  mole  fraction 

z  distance  measured  parallel  to  anode  surface  in  the 

gas  channel  (m) 

Greek  letters 
a  lattice  direction 

0  porosity 

Q  collision  term 


processes  can  be  better  understood  and  thus  optimized  with 
respect  to  the  porous  geometry,  leading  to  improvements  in 
SOFC  performance.  A  detailed  understanding  of  mass  trans¬ 
fer  is  one  of  the  first  steps  in  developing  a  comprehensive 
mathematical  model  of  a  SOFC.  Additional  effects  including 
surface  adsorption/desorption  and  chemical  reactions  (impor¬ 
tant  for  simulating  internal  reforming),  while  important,  are 
beyond  the  scope  of  this  article.  The  two-dimensional  model 
of  the  porous  SOFC  anode  is  shown  in  Fig.  1(b).  The  mole  frac¬ 
tions  of  H2,  H2O  and  N2  are  known  at  the  fuel  or  gas  side  and  the 
mole  flux  of  H2  (and  hence  of  H2O)  is  known  at  the  TPB  based 
on  the  current  drawn  by  the  SOFC.  For  simplicity,  the  TPB  is 
treated  like  an  infinitely  thin  boundary  in  this  study,  although 
the  model  developed  here  can  easily  be  adapted  for  the  case 
where  the  exact  geometric  location  of  the  TPB  is  available  from 
imaging  techniques.  The  solid  obstacles  are  impermeable  to  all 
species  and  no- slip  boundary  conditions  are  applied  at  the  obsta¬ 
cle  surfaces.  Along  the  z  direction,  a  reasonable  height  His  used 
so  that  it  allows  a  realistic  representation  of  the  porous  geometry. 
Periodic  boundary  conditions  are  employed  at  z  =  0  and  z  =  H  for 


2e- 


computational  efficiency.  Because  the  intended  application  is  an 
anode  supported  SOFC,  diffusion  of  O2  through  the  cathode  is 
not  modeled. 

The  lattice  Boltzmann  method  (LBM)  is  used  to  model 
species  transport  through  the  porous  anode.  In  recent  times,  the 
LBM  has  been  well  established  as  an  alternative  computational 
fluid  dynamics  (CFD)  approach  for  a  variety  of  applications 
[1].  The  LBM  has  been  gaining  popularity  owing  to  its  ease  of 
implementation  and  ability  to  capture  the  right  physics  of  flow 
and  diffusion  over  a  wide  range  of  Knudsen  numbers  compared 
to  the  more  traditional  approaches  based  on  the  Navier-Stokes 
equations.  Of  the  many  LBM  variants  designed  to  model  multi- 
component  diffusion,  the  method  chosen  here  is  based  on  the 
two-fluid  model  initially  proposed  by  Luo  and  Girimaji  [2] 
and  extended  by  McCracken  and  Abraham  [3]  in  order  to 
model  components  with  different  molecular  weights.  It  has  been 
shown  that  LBM  can  model  single-component  flow  in  the  non¬ 
continuum  flow  regime  [4,5] .  This  feature  is  important  for  SOFC 
applications  because  high  temperatures  and  small  pore  sizes 
cause  non-continuum  diffusion  effects  to  become  important. 
However,  capability  of  the  LBM  to  model  non-continuum  mass 
diffusion  is  the  focus  of  a  separate  study  within  our  group  and 
this  work  focuses  on  continuum  transport.  The  LBM  is  devel¬ 
oped  for  three  gaseous  species  (1  =  H2,  2  =  H2O,  3  =  N2)  with  no 
phase  changes.  It  is  assumed  that  the  mass  diffusion  occurs  at 
constant  temperature  and  transport  properties  like  binary  diffu¬ 
sion  coefficients  are  constant.  The  LBM  model  used  here  has 
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been  validated  with  the  Stefan-Maxwell  (SM)  equations  for 
simplified  geometries  in  both  ID  and  2D  cases  [6]. 

Modeling  of  mass  transport  in  a  SOFC  has  been  performed 
by  a  large  number  of  research  groups,  [7-15].  However,  in  most 
of  these  prior  studies,  parameters  such  as  the  porosity,  tortuos¬ 
ity,  and  mean  pore  radius  influence  mass  transport  through  the 
porous  medium  and  affect  eventual  predictions  of  cell  perfor¬ 
mance.  While  the  porosity  and  mean  pore  radius  can  be  measured 
experimentally,  tortuosity  is  often  obtained  via  empirical  rela¬ 
tions  that  can  change  from  one  microstructure  to  the  next.  Thus, 
these  models  do  not  have  the  capability  to  model  the  fundamental 
effects  of  electrode  geometry  on  SOFC  mass  transport.  The  goal 
of  the  present  work  is  to  use  the  LBM  approach  to  model  multi- 
component  diffusion  in  complex  geometries  without  relying  on 
empirical  factors  like  tortuosity.  The  relationship  between  mass 
transport  and  concentration  polarization  is  directly  explored  for 
geometries  representative  of  actual  SOFC  anodes. 

The  remaining  part  of  the  paper  is  organized  as  follows.  In 
Section  2,  the  LBM  for  three-component  diffusion  used  in  the 
present  study  is  reviewed.  Parallel  implementation  of  the  LBM 
algorithm  is  discussed  next  with  the  help  of  speed  up  plots. 
Imaging  techniques  used  to  obtain  a  digital  representation  of 
the  anode  geometry  are  described  in  Section  3.  Results  obtained 
from  the  LBM  model,  parametric  studies  and  validation  of  LBM 
model  results  are  discussed  in  Section  4  and  conclusions  are 
presented  in  Section  5. 

2.  Lattice  Boltzmann  method 

The  first  step  in  the  LBM  is  to  overlay  the  solution  domain 
with  a  discrete  lattice  consisting  of  a  network  of  node  points.  A 
uniform,  square  lattice  is  used,  with  obstacles  occupying  part  of 
the  lattice.  The  LBM  essentially  consists  of  two  basic  steps  that 
are  carried  out  at  each  node  point  x  that  is  not  inside  an  obsta¬ 
cle:  streaming  and  collision.  Streaming  represents  movement  of 
particles  of  each  species  i  (where  i=  1,  2  or  3)  along  specified, 
nearest  neighbor  lattice  directions  ela  {a  =  0,  1,  2,  3,  4,  5,  6,  7, 
8}.  These  velocity  directions  are  based  on  the  D2Q9  model  [16]. 
The  collision  term  represents  interactions  between  particles  of 
the  same  or  different  species  as  they  arrive  at  any  given  node. 
These  steps  are  combined  together  in  Eq.  (3),  which  is  called  as 
the  lattice  Boltzmann  equation. 

4(x  +  e^,^+l)-/>,0  =  4  (3) 

In  Eq.  (3),  fla  is  the  velocity  particle  distribution  function 
(PDF)  along  direction  e^,  ns  the  time-level  and  L2la  represents 
the  collision  term  for  species  i  along  the  direction  a.  Physically, 
fla  at  any  given  location  x  represents  the  number  of  particles  of 
species  i  at  that  location  likely  to  be  moving  with  a  velocity  of  . 
Eq.  (3)  essentially  provides  a  time  marching  scheme  to  advance 
the  PDFs  in  time  at  each  lattice  point.  Moments  of  the  PDF  fla 
over  velocity  space  lead  to  the  species  density  and  velocity. 

Boundary  conditions  on  the  PDF  at  the  left  and  right  bound¬ 
ary  of  the  solution  domain  illustrated  in  Fig.  1(b)  are  assigned 
based  on  a  treatment  similar  to  that  used  by  Zou  and  He  [17], 
but  are  modified  such  that  the  average  velocity  before  and  after 


collision  is  used  [18].  Particle  density  (mole  fraction)  is  speci¬ 
fied  at  v  =  0  and  the  x-component  of  particle  velocity  (mole  flux) 
is  specified  at  x-L  for  H2  and  H2O.  For  the  inert  N2,  the  mole 
flux  at  the  TPB  is  zero.  Obstacles  within  the  domain  are  treated 
like  impermeable  solids  with  zero  mass  flux  in  a  direction  nor¬ 
mal  to  their  surface.  The  tangential  velocity  of  all  species  at 
solid  obstacles  is  also  assumed  to  be  zero.  To  implement  these 
conditions  in  the  LBM  model,  bounce-back  conditions  [17] 
are  used  on  lattice  points  corresponding  to  obstacles.  Species 
mole  fractions  are  initially  assigned  the  same  values  over  the 
entire  domain  and  these  are  their  respective  values  at  v  =  0.  The 
v-component  of  species  velocities  (for  H2  and  H2O)  are  ini¬ 
tially  based  on  the  inlet  concentration  and  specified  mole  flux 
at  x-L.  The  y-component  of  velocity  is  initially  zero  every¬ 
where  in  the  solution  domain.  For  lattice  points  that  lie  on  the 
solid  boundary,  both  velocity  components  are  zero.  Using  Eq. 
(3),  the  PDFs  are  then  evolved  in  time  until  a  steady  solution  is 
obtained. 

The  LBM  algorithm  used  in  this  work  [6]  is  implemented 
in  parallel  on  a  SGI  Altix  3700  at  the  University  of  Connecti¬ 
cut.  A  vertical  decomposition  of  the  solution  domain  is  carried 
out  with  each  processor  handling  a  specific  number  of  rows. 
Fig.  2(a)  shows  an  example  anode  geometry  modeled  using  a 
51  x  51  lattice,  with  calculations  split  among  five  processors. 
The  message  passing  interface  (MPI)  library  is  used  for  com¬ 
munication  of  data  between  processors  during  each  time  step. 
In  the  input  geometry  file  shown  in  Fig.  2(a),  solid  obstacles 
are  represented  by  0  and  gas  flow  passages  are  represented  by 
1 .  This  input  file  is  created  using  image  processing  techniques 
discussed  in  Section  3. 

For  a  parallel  simulation,  speed-up  for  N  processors  is  defined 
as  the  ratio  of  run-time  for  a  single  processor  to  the  run-time 
using  N  processors.  Ideally,  if  N  processors  are  used,  the  run¬ 
time  should  fall  by  a  factor  of  N,  leading  to  a  speed-up  of  N.  A 
speed  up  plot  for  typical  LBM  model  runs  is  shown  in  Fig.  2(b). 
The  white  bars  indicate  the  ideal  or  linear  speed-up.  It  can  be 
seen  that  the  actual  speed-up  (indicated  by  the  colored  bars)  falls 
short  of  linear,  especially  for  large  N.  This  occurs  due  to  the 
relatively  large  time  spent  in  communicating  data  between  pro¬ 
cessors.  However,  a  reasonable  speed-up  is  obtained,  allowing 
for  cases  studied  to  achieve  convergence  in  a  reasonable  time. 
It  is  also  quite  encouraging  that  for  large  grid  sizes,  the  LBM 
exhibits  linear  speed  up  through  a  larger  number  of  processors. 
The  current  LBM  code  is  thus  readily  portable  onto  larger  vector 
computers  for  solving  more  complex  problems  using  a  higher 
lattice  resolution. 

3.  Imaging  anode  microstructure 

At  present,  SEM  images  like  the  one  shown  in  Fig.  3  are  first 
used  to  construct  a  detailed  model  image.  This  image  (black 
and  white)  is  similar  to  the  micro  structure  image,  but  has  a  bet¬ 
ter  contrast  between  solid  and  void  spaces.  The  model  image  is 
then  converted  to  a  digital  format  where  solid  (black)  and  void 
(white)  regions  are  represented  by  0  and  1 ,  respectively.  The  con¬ 
version  is  performed  using  ImageJ  (http://rsb.info.nih.gov/ij/). 
In  reality,  the  pores  are  connected  via  pathways  in  and  out  of 
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(b) 


Number  of  processors  (N) 


Fig.  2.  (a)  Vertical  decomposition  of  a  51  x  51  lattice  using  five  processors,  (b)  Speed  up  behavior  of  the  LBM  model  using  vertical  decomposition. 


the  plane  of  the  image,  but  this  information  can  not  be  obtained 
without  using  a  3D  imaging  technique.  Thus,  while  construct¬ 
ing  the  model  image,  it  is  sometimes  necessary  to  artificially 
change  the  digitized  image  to  provide  a  connected  network  of 
pores  and  thus  have  clear  pathways  for  mass  diffusion.  In  addi¬ 
tion  to  the  above  limitation,  the  anode  materials  used  for  SEM 
(obtained  from  actual  SOFC  prototypes)  are  of  a  fixed  porosity 
value  (~30%)  and  geometry  data  for  different  porosities  is  not 
available.  In  order  to  carry  out  parametric  studies,  the  porosity  is 
reduced  by  enlarging  the  size  of  solid  obstacles  (reducing  void 
spaces)  and  by  adding  more  obstacles  in  the  domain  if  neces¬ 
sary.  Because  of  these  limitations,  the  2D  geometry  should  be 
viewed  as  qualitatively  resembling  actual  SOFC  anodes.  A  more 
accurate  representation  will  have  to  wait  until  more  advanced 
3D  imaging  techniques  are  available.  The  lattice  size  for  the 
FBM  will  need  to  be  large  enough  to  adequately  resolve  the 
microstructure.  At  the  same  time,  it  should  lead  to  reasonable 
run-times.  Run  time  foral51  x  151  lattice  using  25  processors 
was  about  135  min. 


4.  Results  and  discussion 


The  two  main  dimensionless  parameters  governing  contin¬ 
uum  mass  transport  through  a  porous  medium  are  the  porosity 
0,  which  is  defined  as  the  fraction  of  anode  area  (volume  in  3D) 
occupied  by  void  space,  and  the  dimensionless  flux  (J*) 


J* 


JL 

Ct£>3i 


(4) 


In  Eq.  (4),  J  is  the  mole  flux  at  the  TPB,  L  is  the  length  of 
the  medium,  Op  is  the  total  molar  concentration,  and  D31  is  the 
binary  diffusion  coefficient  between  H2  and  N2.  Physically,  / * 
represents  the  relative  strength  of  convective  and  diffusive  mass 
transport  of  hydrogen  gas  through  N2.  Because  mass  transport  in 
a  SOFC  tends  to  be  diffusion  dominated,  values  of  f  tend  to  be 
quite  small,  typically  of  the  order  of  10-2.  FBM  parameters  can 
be  arbitrarily  chosen  as  long  as  the  corresponding  dimension¬ 
less  numbers  are  correct.  For  example,  to  simulate  a  case  where 
f  =  0. 16,  the  FBM  parameters  (in  lattice  units)  can  be:  L  =  150, 
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Digital  image  file 
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Model  image 


Fig.  3.  Construction  of  digital  geometry  input  file  from  a  SEM  image. 


J  =  0.069,  CT  =  600,  Dn  =  0.337, £>23  =  0.069  and  D31  =  0. 1085. 
Unless  indicated  otherwise,  mole  fractions  of  H2,  H20  and  N2 
are  0.47,  0.03  and  0.5,  respectively  at  v  =  0.  Because  it  involves 
the  mole  flux,  J*  can  be  related  to  the  current  density  of  the  fuel 
cell  via  [19] 


i  =  2F  x  J 


(5) 


In  Eq.  (5),  i  is  the  current  density  and  F  is  Faraday’s  con¬ 
stant.  In  the  present  study,  the  mole  flux  remains  constant 
over  the  TPB,  leading  to  a  constant  current  density  at  the 
anode-electrolyte  interface.  The  current  is  non-dimensionalized 
using  f  =  iLKflFCjD^i )  =  f .  The  voltage  or  EMF  developed  by 
the  cell  E  can  be  obtained  using  the  fuel  (H2),  product  (H20) 
and  oxygen  (02)  partial  pressures  (or  mole  fractions)  at  the  TPB 
in  the  Nernst  equation  [19] 


n  RT 

E°  + -  In 

2  F 


Pn2o 


(6) 


In  Eq.  (6),  is  the  EMF  at  standard  pressure  and  R  is  the  gas 
constant.  As  mentioned,  diffusion  through  the  cathode  is  not 
modeled  and  the  partial  pressure  of  02  at  the  cathode  TPB  is 
assumed  to  be  a  constant  21%.  The  concentration  polarization 
(A  Vconc)  is  defined  to  be  the  loss  of  voltage  due  to  the  drop  in 
fuel  (H2)  concentration  across  the  anode  from  the  gas  channel 
v  =  0  to  the  TPB  x  =  L. 


A  Vconc  —  ^gas  channel  Etpb  (7) 

In  Eq.  (7),  £gas  channel  is  evaluated  using  H2  and  H20  mole  frac¬ 
tions  of  0.47  and  0.03  respectively,  while  Etpb  is  evaluated 


based  on  the  average  H2  and  H20  mole  fractions  obtained  at 
x-L  from  the  LBM  model. 

4. 1.  Validation  of  the  LBM  model 

Experimental  data  for  just  concentration  polarization  is 
scarce,  so  the  LBM  results  have  been  compared  with  some 
prior  models  that  calculate  concentrcaation  polarization  (Chan 
et  al.  [8],  Zhao  and  Virkar  [10]).  In  order  to  have  a  fair  com¬ 
parison,  the  dimensional  current  i  used  in  these  prior  models 
needs  to  be  converted  to  a  non-dimensional  current  density 
i* .  This  conversion  is  illustrated  for  one  particular  case.  The 
following  data  for  the  anode  is  obtained  from  Chan  et  al. 
[8]:  porosity  =  30%,  tortuosity  =  6 ,  L  =  750  x  10-6  m.  Based  on 
Eq.  (5),  a  current  density  i  =  20,000 Am-2  corresponds  to  a 
mole  flux  J- 0.1036  mol  m-2  s-1.  The  binary  diffusivity  is 
£>13  =  1.085  x  10_4m2  s-1  andCx  =  11.4  mol  m-3  at  a  pressure 
P=  1  atm  and  temperature  T=  1073  K  based  on  the  ideal  gas  law. 
Note  that  tortuosity  is  not  a  parameter  in  the  LBM  model,  but 
is  an  empirical  parameter  used  to  change  the  binary  diffusivity 
value  by  Chan  et  al.  [8],  Zhao  and  Virkar  [10]  and  by  many 
other  models  that  are  based  on  a  ID  representation  of  the  porous 
medium.  The  above  set  of  parameters  {L,  /,  Cr,  £>i3}  lead  to 
/*  =  f  =  0.0628.  Thus,  a  current  density  i  =  20,000  Am-2  in  the 
study  by  Chan  et  al.  [8]  corresponds  to  a  non-dimensional  cur¬ 
rent  density  i*  =  0.0628.  This  relationship  is  used  to  convert  the 
data  provided  by  Chan  et  al.  [8]  (Fig.  3  from  their  paper)  to 
dimensionless  form.  A  similar  conversion  is  performed  on  the 
data  provided  by  Zhao  and  Virkar  [10]  in  Fig.  21  of  their  paper. 
This  latter  study  has  a  porosity  of  32%,  which  is  fairly  close  to 
the  30%  porosity  data  calculated  using  the  LBM  model. 
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Fig.  4.  Validation  of  the  LBM  model  for  a  porosity  of  30%. 


Fig.  4  compares  the  concentration  polarization  from  the  LB 
model  with  data  from  Chan  et  al.  [8]  and  Zhao  and  Virkar  [10]. 
The  LBM  polarization  curve  corresponds  to  the  30%  porosity 
data  from  Fig.  8,  with  many  additional  data  points.  It  can  be 
observed  that  Chan  et  al.  [8]  predict  an  initially  convex  and  ulti¬ 
mately  concave  shape  of  the  concentration  polarization  curve. 
However,  the  limiting  current  density  in  the  LBM  model  is 
lower,  so  the  entire  portion  of  the  concave  part  of  the  curve 
(for  large  f)  is  not  obtained  at  30%  porosity.  This  concave 
region  is  also  absent  in  the  results  of  Zhao  and  Virkar  [10].  It 
is  observed  the  LBM  model  predictions  are  fairly  close  to  these 
two  prior  models,  indicating  that  the  LBM  approach  is  accurate 
enough  to  model  concentration  polarization  in  a  SOFC  using  a 
2D  model. 

The  differences  observed  in  the  predicted  concentration 
polarization  may  be  because  of  the  2D  geometry  used  in  this 
work,  non-continuum  effects  and/or  because  of  empirical  fac¬ 
tors  like  tortuosity  used  in  prior  models.  The  electrochemical 
coupling  between  current  density  and  species  concentration  may 
be  important  when  the  limiting  current  condition  is  reached.  In 
addition,  it  should  also  be  pointed  out  that  the  LBM  model  pre¬ 
sented  here  uses  a  ternary  system  {H2,  H2O,  N2},  while  both 
Chan  et  al.  [8]  and  Zhao  and  Virkar  [10]  employ  a  binary  system 
of  H2  and  H2O. 
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Fig.  5.  Mole  fraction  distribution  for  H2  (a),  H2O  (b)  and  N2  (c)  inside  a  porous 
anode  for  0  =  0.5  and  J*  =0.16. 


4.2.  Mass  transport  through  porous  media 

Fig.  5  illustrates  a  typical,  steady  state,  mole  fraction  varia¬ 
tion  of  H2,  H2O  and  N2  inside  the  porous  anode  for  0  =  0.5  and 
J*  =  0.16.  A  high  value  of  J*  =0.16  is  used  in  order  to  clearly 
indicate  mole  fraction  variation  of  H2,  H2O  and  N2  in  a  typical 
porous  geometry.  Smaller  values  of  J*  do  not  produce  apprecia¬ 
ble  changes  in  species  mole  fractions.  In  general,  diffusive  mass 
transport  of  any  species  occurs  from  higher  to  lower  concentra¬ 
tions.  Because  H2  is  consumed  at  the  TPB,  the  mole  fraction  of 
H2  at  the  TPB  is  lower  compared  to  the  gas  side.  Similarly,  as 
H2O  is  generated  at  the  TPB,  mole  fraction  of  H2O  at  the  TPB 
is  higher.  For  the  inert  N2,  even  if  there  is  no  mass  transfer,  the 
mole  fraction  at  the  TPB  reduces.  For  all  three  species,  mole 
fraction  at  the  TPB  varies  with  z  depending  on  the  local  struc¬ 
ture  of  the  porous  medium.  Therefore,  mole  fractions  of  H2  and 
H2O  at  the  TPB  are  averaged  along  z  for  parametric  studies  and 
for  calculating  the  concentration  polarization. 


4.3.  Effect  of  porosity  and  dimensionless  flux  on  mass 
transport 

The  porosity  in  the  LB  model  is  controlled  by  changing 
the  amount  of  void  spaces  present  in  the  solid  matrix.  Anodes 
of  different  porosities  from  0  =  0.25-0.5,  constructed  using  a 
151  x  151  lattice,  are  shown  in  Fig.  6.  For  a  given  J* ,  as  the 
porosity  reduces,  there  are  less  diffusion  pathways  through 
the  domain  and  a  larger  concentration  difference  is  required 
in  order  to  maintain  the  same  amount  of  molar  flux.  As  the 
mole  fraction  of  H2  is  fixed  at  x  =  0,  the  mole  fraction  of  H2 
at  the  TPB  reduces  as  the  porosity  reduces.  The  mole  fraction 
values  at  a  porosity  of  100%  correspond  to  the  case  where 
there  are  no  obstacles.  As  a  check,  it  is  verified  that  predicted 
values  using  the  LB  model  match  those  obtained  from  the  ID 
solution  to  the  SM  equations.  For  a  fixed  porosity  value,  as  J* 
is  increased,  concentration  gradient  increases  to  accommodate 
the  larger  mole  flux.  Due  to  increased  current  draw,  H2  mole 
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Fig.  6.  Anode  geometries  used  to  study  the  effect  of  porosity  (f)  with  a  151  x  151  lattice. 


fraction  at  the  TPB  reduces.  If  /  is  too  large,  the  limiting 
current  condition  is  reached  because  the  H2  concentration  at 
the  TPB  tends  to  approach  zero.  This  effect  is  magnified  at 
low  porosity  values.  Thus,  a  limit  is  imposed  on  the  magnitude 
of/. 

Both  of  these  effects  are  illustrated  in  Fig.  7,  which  plots  the 
average  H2  mole  fraction  at  the  TPB  as  a  function  of  poros¬ 
ity  for  different  /  values.  It  can  be  seen  that  for  high  /,  the 
porosity  has  to  be  above  a  threshold  value  for  solutions  to 
be  obtained.  The  variation  with  porosity  is  highly  non-linear, 
except  for  the  trivial  case  of  /  =  0,  representing  no  current 
draw.  This  study  assumes  that  the  current  density  at  the  TPB 
is  independent  of  species  concentration  at  that  point.  A  detailed 
analysis  is  currently  in  progress  in  order  to  establish  the  nature  of 
the  interdependence  of  the  diffusion  and  electrochemical  reac¬ 
tion  processes.  Once  this  analysis  is  complete,  the  LBM  model 
can  be  suitably  adapted  to  take  this  coupling  into  account  if 
required. 


4.4.  Effect  of  porosity  and  dimensionless  flux  on 
concentration  polarization 

Once  the  average  H2  and  H2O  mole  fractions  at  the  anode 
TPB  are  known  from  the  LBM,  anode  concentration  polarization 
curves  can  be  constructed  using  Eq.  (7)  to  study  SOFC  perfor¬ 
mance.  Because  lower  fuel  concentration  at  the  TPB  leads  to 
higher  polarization  losses,  it  can  be  deduced  that  these  losses 
are  more  for  low  porosity  values  and  high  /  (larger  current 
density).  Concentration  polarization  curves  are  calculated  for  a 
range  of  i*  values  and  for  different  anode  porosities  in  Fig.  8.  It 
is  found  that  the  polarization  curve  is  almost  linear  at  interme¬ 
diate  f\  with  steeper  gradients  at  the  lower  and  higher  ends.  As 
expected,  since  f  is  proportional  to  /,  the  lowest  polarization 
occurs  at  small  1  . 

The  general  trend  of  the  polarization  curve  is  similar  for  dif¬ 
ferent  anode  porosities.  However,  for  a  given  current  density 
f\  as  the  porosity  is  reduced,  the  polarization  increases.  This 
is  due  to  the  reduction  in  H2  mole  fraction  at  the  anode  TPB 


J*  =  0 

-m-  J*  =  0.032 
—A—  J*  =  0.064 
-D-r  =  0.096 
=  0.128 
J*  =  0.160 
J*  =  0.192 


Fig.  7.  Effect  of  porosity  and  J*  on  H2  mole  fraction  at  the  TPB .  LBM  parameters 
(in  lattice  units):  L-  150,  J- 0.069,  D\^  =  0.337,  £>23  =0.069  and  D31  =0.1085. 
Mole  fractions  of  H2,  H2O  and  N2  are  0.47,  0.03  and  0.5  respectively  at  x  =  0. 
f  changed  by  changing  the  value  of  Cj. 


Non-dimensional  current  density  ( i* ) 

Fig.  8.  Anode  concentration  polarization  for  different  porosity  values. 
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at  low  porosity.  The  maximum  current  drawn  by  the  SOFC  is 
also  smaller  for  a  low  porosity.  In  Fig.  8,  the  increasing  slopes 
of  the  polarization  curves  indicate  that  the  polarization  is  more 
sensitive  to  changes  in  current  density  at  low  porosities.  Some 
curves  in  Fig.  8  stop  at  a  lower  voltage  because  solutions  are 
not  obtained  at  the  next  i*  value  due  to  the  H2  mole  fraction 
becoming  zero  (limiting  current  condition). 

Increasing  the  anode  porosity  leads  to  better  mass  transport 
of  fuel  and  lower  concentration  polarization.  However,  a  high 
porosity  can  reduce  the  number  of  electrochemical  reaction  sites 
as  well  as  the  pathways  for  the  flow  of  electrons  through  the 
anode,  thereby  increasing  activation  and  ohmic  losses  (not  mod¬ 
eled  here).  High  porosity  can  also  lead  to  reduced  mechanical 
strength.  An  optimum  porosity  and  micro  structure  geometry 
may  exist  where  the  combined  polarization  losses  (concen¬ 
tration  +  activation  +  ohmic)  are  a  minimum.  In  practice,  the 
porosity  is  often  comparatively  larger  near  the  gas  channel  and 
is  lower  near  the  electrolyte  to  provide  more  surface  area  for 
electrochemical  reactions  to  take  place.  Using  the  LBM  model 
developed  here,  SOFC  anode  microstructure  can  be  optimized 
to  give  the  best  performance. 

5.  Conclusions 

The  lattice  Boltzmann  method  (LBM)  is  used  as  a  tool  to 
model  mass  transport  of  H2  (fuel)  and  H2O  (product)  in  the 
presence  of  N2  (inert  gas)  in  the  porous  anode  structure  of  a 
SOFC.  The  porous  structure  used  in  the  LBM  model  is  based  on 
SEM  images,  which  are  converted  to  digital  form.  Results  for 
mole  fraction  distributions  are  obtained  for  a  range  of  porosity 
values  and  dimensionless  flux  values.  For  high  porosity,  predic¬ 
tions  using  the  LB  model  approach  those  obtained  from  solving 
the  continuum,  ID  Stefan-Max  well  equations.  As  porosity  is 
reduced,  the  fuel  concentration  at  the  TPB  rapidly  falls  and  even¬ 
tually  approaches  zero.  For  any  porosity,  there  is  a  condition  of 
maximum  J*  above  which  a  limiting  current  density  in  the  SOFC 
can  be  predicted. 

Advantages  of  the  LBM  model  are  that  a  detailed  analysis  of 
mass  transfer  can  be  carried  out  for  an  actual  anode  microstruc¬ 
ture  and  fuel  cell  performance  can  be  judged  by  calculating 
the  concentration  polarization.  Complex  geometries  can  be  eas¬ 
ily  handled  and  both  continuum  and  non-continuum  diffusion 
through  the  pores  can  be  modeled.  Because  of  the  level  of  detail, 
effect  of  anode  geometry  on  voltage  polarization  can  be  directly 
studied  and  the  geometry  optimized  to  give  the  best  performance. 
Unlike  prior  models,  empirical  parameters  like  tortuosity  are  not 
used  as  fitting  parameters  to  account  for  the  porous  medium. 

A  limitation  of  the  present  study  is  the  absence  of  three- 
dimensional  (3D)  connectivity  of  the  anode  pore  structures,  but 


this  will  be  overcome  in  the  future  as  better  image  process¬ 
ing  techniques  start  becoming  available.  Some  recent  image 
processing  work  [20]  is  a  step  in  the  right  direction,  but  non¬ 
destructive  imaging  techniques  like  XMT  that  will  minimize 
disruption  of  the  anode  micro  structure  also  need  to  be  devel¬ 
oped  further.  The  LBM  model  itself  can  be  adapted  to  3D  and 
boundary  conditions  can  be  changed  such  that  species  genera¬ 
tion  and  consumption  can  be  applied  at  specific  TPB  locations 
within  the  microstructure.  In  addition  to  mass  transport,  acti¬ 
vation  and  ohmic  losses  in  the  anode  will  be  considered  in 
future  studies  as  3D  imaging  data  becomes  available.  Incor¬ 
poration  of  these  effects  in  the  present  LBM  model  will  lead 
to  a  highly  accurate  and  useful  tool  for  SOFC  design  and 
analysis. 
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